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ABSTRACT 

We study the self-consistent, linear response of a galactic disc to non-axisymmetric 
perturbations in the vertical direction as due to a tidal encounter, and show that the 
density distribution near the disc mid-plane has a strong impact on the radius beyond 
which distortions like warps develop. The self-gravity of the disc resists distortion 
in the inner parts. Applying this approach to a galactic disc with an exponential 
vertical profile, Saha & Jog showed that warps develop beyond 4-6 disc scalelengths, 
which could hence be only seen in HI. The real galactic discs, however, have less 
steep vertical density distributions that lie between a sech and an exponential profile. 
Here we calculate the disc response for such a general sec/i 2 /™ density distribution, 
and show that the warps develop from a smaller radius of 2-4 disc scalelengths. This 
naturally explains why most galaxies show stellar warps that start within the optical 
radius. Thus a qualitatively different picture of ubiquitous optical warps emerges for 
the observed less steep density profiles. The surprisingly strong dependence on the 
density profile is due to the fact that the disc self-gravity depends crucially on its 
mass distribution close to the mid-plane. General results for the radius of onset of 
warps, obtained as a function of the disc scalelength and the vertical scaleheight, are 
presented as contour plots which can be applied to any galaxy. 

Key words: Galaxies: kinematics and dynamics - Galaxies: spiral - Galaxies: struc- 
ture 



1 INTRODUCTION 

It is a common knowledge now that the galaxies are not 
isolated structures in the universe, rather galaxy interac- 
tions including mergers are common. Hence understanding 
the dynamics of these interactions and the effects they pro- 
duce on the structure and evolution within galaxies has been 
of much interest. One such effect is the generation of asym- 
metric features due to tidal encounters between galaxies. 

Spiral galaxies are known to display a variety of non- 
axisymmetric features, both in the plane and also in the 
direction perpendicular to the plane. The most common ver- 
tical distortion is a warp, a feature of azimuthal wavenumber 
m = 1. Most nearby galaxies show such integral-sign or s- 
shaped warps in their outer parts (e.g., Binney & Tremainc 
1987). A similar planar distortion commonly seen is the lop- 
sidedness in disc corresponding to azimuthal wavenumber 
m — 1 (Jog & Combes 2009). The origin of warp is not yet 
clear despite a long search. A commonly suggested mech- 
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anism for the origin of warps is due to the tidal interac- 
tion with its neighbours (e.g., Schwartz 1985, Zaritsky & 
Rix 1997). Weinberg (1995) studied the generation of warp 
in our Galaxy due to perturbation from the neighbouring 
Large Magellanic Cloud. 

Vertical distortions other than warps are also commonly 
observed in spiral galaxies. Small-scale corrugations in exter- 
nal galaxies have been studied (Quiroga 1984) and so have 
been the distortions in the stellar distributions of galaxies 
(Florido et. al 1991). Scalloping in HI in the outer regions of 
our Galaxy has been studied too (Kulkarni, Blitz & Heiles 
1982). Recently Matthews & Uson (2008) find that corruga- 
tion in HI is seen in IC2233 even in the inner regions. 

While warps are seen mostly in the outer parts of a 
galactic disc, surprisingly little work has been done to dis- 
cuss the radius at which warps develop. Saha & Jog (2006) 
proposed this as being determined due to the self-consistent 
disc response to a tidal field. The disc self-gravity resists dis- 
tortion in the inner parts, and only in the outer parts does 
the disc begin to respond to the external potential. 

In this paper, we continue with this approach and study 
self-consistent vertical distortions for different forms of ver- 
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tical density distributions in the disc. We are motivated by 
the fact that in the previous studies, not much attention has 
been paid to the effect the form of density distribution for 
the responding disc might have on the overall behaviour of 
the system. For the purpose of modeling, the vertical distri- 
bution was taken to be exponential for mathematical sim- 
plicity in an earlier work (Saha & Jog 2006). 

The study of vertical distribution of stars in galactic 
discs has an interesting history. The vertical distribution 
was at first deduced to be of type seek 2 as resulting for 
an isothermal disc (Spitzer 1942). However, observations 
showed a steeper profile closer to a seek or an exponen- 
tial both for our Galaxy (e.g., Gilmore & Reid 1983, Kent 
et al. 1991) as well as for external galaxies (e.g., Wainscoat, 
Freeman & Hyland 1989, Rice et al. 1996). An exponen- 
tial profile all the way to the mid-plane is unphysical and 
hence van der Kruit (1988) proposed a generalized function 
sech 2 /™ to represent the vertical density distribution. In this 
scheme, n = 1 and 2 correspond to a sech 2 and a sech dis- 
tribution and as n tends to oo it asymptotically approaches 
an exponential distribution. Later observers analyzed their 
data and cast in this format, and have shown that a true 
density distribution is less steep than an exponential and in 
most case lies between a sech and an exponential distribu- 
tion (Barteldrees & Dettmar 1994, de Grijs, Pcletier & van 
der Kruit 1997). While dust extinction prevents a determi- 
nation of the stellar density profile close to the mid-plane, 
the near-infrared bands do not have this limitation and rep- 
resent the true density profiles representing the old stars. 
For a sample of 24 galaxies studied in the K s band, de Grijs 
et al. (1997) find that, a mean value of < 2/n >= 0.5 cor- 
responding to the n index = 4. Thus the vertical density 
profile for stars in a typical galactic disc lies between a sech 
and an exponential profile, being closer to a sech. 

In this paper, we study the vertical response of an ax- 
isymmetric disc to an external imposed perturbation, where 
the disc density follows such a general sech 2 ^ n distribution. 
We also study the response of an exponential disc for the 
sake of comparison. The motivation for our work comes from 
the fact that the matter distribution close to the mid-plane 
contributes strongly to the vertical self-gravitational force 
(Banerjee & Jog 2007), hence it is plausible that the differ- 
ent vertical density profiles affect the disc response in dif- 
ferent ways. We show that indeed the disc response has a 
strong dependence on the vertical density distribution. This 
in turn significantly affects the radius for the onset of various 
non-axisymmetric features along the vertical direction. 

Section 2 of the paper presents the details of the model 
and the methods of calculation, while section 3 presents the 
results. Section 4 presents the discussion and section 5 con- 
cludes the paper. 



2 DETAILS OF THE MODEL 

2.1 Density response of the perturbed disc 

In this paper, we study the linear response of the disc to 
an external imposed perturbation potential for a general 
sechP' 71 ' vertical distribution, as well as the exponential 
case for comparison (see Section 1). In an earlier work, 
the density distribution perpendicular to the plane of the 



disc was taken to be exponential for simplicity (Saha & Jog 
2006). Here we build on the formulation developed in that 
earlier paper, with the main difference being that the den- 
sity profile p(z) is taken to be a general, flat-core of type 
sech ( - 2/n) in this case. Also, the radial variation of the po- 
tential is taken into account properly here, which leads to a 
more realistic behaviour of the response function. 

We use cylindrical co-ordinates. The stellar density dis- 
tribution of the unperturbed axisymmetric disc is taken as: 



p d (R,z) = h{R)h{z) 



(1) 



where fi(R) is a function corresponding to the radial de- 
pendence and fa (z) is a function which takes care of the de- 
pendence in the vertical direction. The radial dependence is 
taken to be exponential (Freeman 1970), while a sech^ 2 ^ n ^ 
dependence in the vertical direction is taken, as based on 
the work by van der Kruit (1988). However, we want to con- 
sider the case when the mid-pane density at a given radius 
is same for any n value, hence we use the form as in Kuijken 
& Gilmore (1989, see eq.[A13]) : 

h{R)=p do exp { - R/RD) ; f 2 (z) = sech 2/ "(^) (2) 

where pdo is the central, mid-plane density and Rd is the 
exponential disc scalelength. Note that for n = 1,2, the func- 
tion /2(z) corresponds to a sech 2 and a sech profile respec- 
tively, while for n— > oo, it reduces to a form exp (—z/zo). 
Thus zo is taken to denote the vertical scaleheight. Hence, 



(3) 



p d (R,z) = p d0 ezp<-™ S ech 2/n (^) 

2Z0 

In order to see the effect that the halo dark matter distri- 
bution might have on the response, for the sake of com- 
pleteness, we also include an axisymmetric spheroidal sys- 
tem with a pseudo-isothermal density profile (de Zeeuw & 
Pfenniger 1988): 



Ph(R,z) 



PhO 



1 + [R? + (zyq 2 )]/R 2 



(4) 



where p h0 , R c and q are the central density, the core ra- 
dius, and the flattening parameter or the vertical to planar 
axes ratio respectively. The halo symmetry plane (z=0) is 
assumed to coincide with the disc plane. 

The density distributions of the halo and the disc are 
related to the corresponding potentials by the Poisson equa- 
tion 



R dR dR 



■ AivGp 



(5) 



where G is the gravitational constant and the total density 
p is: 



p = Pd + Ph 



(6) 



For a "disc plus halo" system, the total potential is the 
combination of the individual potentials of the disc and halo 



(7) 



For a vertical distribution as in eq. (2), the correspond- 
ing integral obtained after inverting the Poisson equation is 
not solvable analytically. However, it can still be represented 
as a sum of integrals which converge fairly rapidly: 



V d {R,z) = -2tvG / dk.J (kR)p d a(a 2 + k 2 ) 
Jo 



-3/22l+(2/n) 
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E 



-2/n\ (1/foXl + wm)e- fc|z| - ke^ +nm)z/z ° 
(l + nm) 2 (l/z ) 2 -k 2 



(8) 



where a — 1/Rd, and Jo(kR) is the cylindrical Bessel func- 
tion of the first kind, of order zero. The m th binomial co- 
efficient of power -(2/n) is denoted by (^ /n ) (Sackett & 
Sparke 1990). If it so happens that n, k and zo have val- 
ues such that the denominator has a zero value, such terms 
can still be evaluated by removing the singularity using 
L' hospital's Rule. The resulting term is then 2 -^— (~^ n ) 
(l + fc|z|)e- fc|z|(2/n) . 

The potential corresponding to the halo density distri- 
bution is not analytically obtainable either, but following 
the more general expression for obtaining potential corre- 
sponding to any distribution, we can reduce it down to a 
form which is numerically integrable (Binney & Tremaine 
1987, eq. 2.88). We have, then 



Vh{R,z) = 2-KG Ph oqRc 

-1/9 2 



R 



Ml + j^i^T^ + z)]dx/{e'x< + 1) 



(9) 

where e = (1 — q 2 ) 1 ^ 2 is the eccentricity of the halo. 

We consider the zero-forcing frequency perturbation po- 
tential due to a perturber situated at a distance D from the 
centre and at an angle of inclination i from the disc mid- 
plane to be: 



Vi(R, (j>, z) = * p0 (i?)(l + — sin i) cos (mq 



(10) 



In the limit that R/D < 1 and z/D 1, * p0 = GM P R/D 2 , 
is a slowly varying function of the radius. We treat this more 
realistic case with the effect of perturbation increasing with 
radius. In the earlier study (Saha & Jog 2006) this was taken 
to be a constant. 

The total potential that a particle in the disc experi- 
ences is a superposition of the potential generated by the 
disc plus the dark matter halo system and the external per- 
turber, given by equations (7) and (10) respectively: 



* total = *(R,z) + V l (R,<l>,z) 



(11) 



To obtain the vertical response of a disc to an imposed per- 
turbation, one needs to solve the equation of motion as a 
forced oscillator, and combine with the continuity equation, 
and with the Poisson equation for the perturbation. See Saha 
& Jog (2006), section 2.1 for details. This gives the following 
for the perturbed motion of a particle: 



H(R)cos{mQ,t) 



where 



H(R) 



8*1 i 



(12) 



(13) 



Here vq and f2o are the vertical and planar frequencies of 
oscillation in the unperturbed system and are defined re- 
spectively as: 
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The perturbed mass density pi is given by: 



pi(R,cj>,z) = -^f-H(R )cos(mq 
oz 



(15) 



where p is the net unperturbed density = Pd + Ph (eg. [6]). 
Using the equation (2) for the vertical density distribution in 
the disc, and noting that dpd/dz = —[pd tanh(nz/2zo)] / z () , 
the perturbed density reduces to: 



Pi(R,<t>,z) = [p d (R,z)tanh(nz/2z ) + Ch(R,h)][ 



zo 

xcos(mq 



(16) 



On simplifying H(R) using eqs. (10) and (13), this reduces 
to: 



pi (R, 4>, z) = p r (R, z) cos(m(f>) 



(17) 



where p r (R,z), the magnitude of the perturbed density, is: 
(GM P R/D 3 )sini 



pr(R, Z 



vo 2 - ^o 2 
[p d (R, z)tanh(nz/2z ) + (h(R, h)} 



(18) 



Here (,h{R,h) (eq.[16]) is defined to be = —zodph/dz, note 
that this is a positive definite quantity for any realistic, cen- 
trally concentrated gravitating system. Further, since the 
perturbation potential has a positive gradient (d^i/dz > 
0), and vq 2 the vertical frequency square is greater than Qo 2 , 
the planar frequency square in an oblate geometry, hence 
the function H(R) > 0. Thus the perturbed disc response 
follows the perturbing potential. This has important impli- 
cations since the corresponding response potential opposes 
the imposed potential, as discussed in the next section. This 
concept was first noted and discussed for the planar pertur- 
bation such as lopsidedness, by Jog (1999). 

Most studies of galactic structure assume a spherical or 
an oblate halo for which the above relation {y{ 2 > fto 2 ) is 
valid (see e.g., the model for our Galaxy by Mera, Chabrier 
& Schaeffer 1998). A prolate halo would not permit this, 
especially in the outer parts where the halo dominates but 
such a halo is not indicated. Recent studies of halo shape ob- 
tained by modeling the observed HI scaleheights have shown 
that the halo is spherical as in our Galaxy (Narayan, Saha 
& Jog 2005, Kalberla et al. 2007), and in UGC 7321 (Baner- 
jee, Matthews, & Jog 2010), or is oblate as in M31 (Banerjee 
& Jog 2008). Here we study a nearly spherical halo, with a 
small oblateness (q, the axis ratio = 0.95). 



2.2 Disc response potential 

We next obtain the potential corresponding to the density 
response of the disc to the imposed potential (eq.[16]). This 
is obtained by inverting the Poisson equation. The Poisson 
equation is given by: 



V 2 ^res P (R,<P,z) 



4TvG Pl (R,(t>,z) 



(19) 



For a disc of non-zero thickness, we solve the Poisson 
equation by using the Green's function approach (Jackson 
1975), see Saha & Jog (2006) for details. This is an integral 
function approach which gives a numerically tractable form 
of the resulting response potential. The Green's function in a 
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cylindrical form for a finite thickness is developed following 
the treatment in Binney & Tremaine (1987, see chapter 2). 

Vres P (R,<t>,z) = -G ( (20) 

JV\r- r'\ 

where V denotes the volume of the galaxy, and 

\r - r\ = [R 2 + R' 2 - 2RR'cos{(t> 1 - <j>) + (z - z') 2 ] 1 ' 2 (21) 

Using the above, the equation for response potential for the 
Fourier component m becomes 



Iv i R2 + 



p r (R', z')cos(m<f>')R' dR'd(j>'dz' 



v [R 2 + R' 2 - 2RR' cos(<j}' -<t>) + (z- z') 2 ] 1 



/2 



(22) 



where p r (R, z) is the magnitude of the density response (see 
eq.[18]). The denominator in the above equation can be sim- 
plified as: 

[R 2 + R' 2 - 2RR' cos(^' - <f>) + (z - z') 2 ] 1/2 
= [(R + R') 2 + (z - z') 2 ] 1 ' 2 x [1 - fcW(m[0 - 0']/2)] 1/2 

(23) 

where fc 2 is defined as 

k 2 = (24 ) 

(R + R') 2 + (z-z') 2 + z s 2 { ' 

Here z s is the small softening parameter added to prevent 
the response potential from blowing up at z — z . 

Next, eq.(22) is simplified by substituting from eq.[23] 
into it: 

p r (R,z)R'dR'dz' 



OO f oo 



<fTes P {RA,z) = -G 



ri 

i 



„ {(R + R') 2 + (z-z>) 2 ]i/ 2 

cosiruf)' d(f>' 
[1 - k 2 cos 2 (m[</> - <t>']/2)} 1 / 2 



(25) 



In the above equation, the integral over (f> can be solved by 
introducing a variable /3 defined as 



p={4>-cj>')/2 +7T/2 

The integral over <f> then reduces to: 

cosrruf)' d(f>' 
[1 - k 2 cos 2 (m[(j> - 0']/2)] 1 / 2 



(26) 



= ±4:Cosm(p x 



I COi 

Jo ~ 



cos2m/3dl3 



k 2 sin 2 p 



Qodd(m,k) 



J 

Jo 



T/2 sin 2 (mj3) 

\l-k 2 sin 2 {P)Y/ 2 P 



(29) 



and K(k) is the complete elliptical integral of the first kind 
defined as: 



K(k) 



tt/2 



[l-Psin 2 ^)] 1 / 



(30) 



where k is as defined in eq.(24). 

In case of even m values, the equation retains a similar 
form except we have, instead of o dd, a similar function 
Oeven defined as 



'even 



r/2 



cos 2 (rn/3) 



{l-k 2 sin 2 pyi 



(31) 



This in turn gives us the equation for the disc response po- 
tential for the even m values as: 



^7elp{R,<i>,z,m) = +4Gcos(m<? 



-,even 
esp 
POO p oo 



p r (R',z')R'dR'dz' 



[{R + R') 2 + 



( g /_z)2]l/2 [2e e ,en(m,fc) - R 



In equations (28) and (32), the quantity p r is equal to 
the magnitude of the response density as given by eq. [18]. 

From here on, our treatment differs from that of Saha 
& Jog (2006), in that we do not take the perturbation po- 
tential, \I/po = GMpR/ D 2 outside the integral over R (which 
was not justified since R ~ R'); and we do not ignore the 
contribution of the halo term within the integral over radius. 
These mathematical refinements, in particular the correct 
treatment of the radial variation in the numerical calcula- 
tions, lead to substantially different results at large radii and 
for high m values compared to Saha & Jog (2006) even for 
the exponential case treated in that paper, see Section 4. 

We obtain the response potential numerically, and find 
that it has a sign opposite to the imposed potential (eq.[10]), 
thus it is negative. Thus the potential corresponding to the 
disc response to the imposed potential alone opposes it. This 
is a general result for any self-gravitating system, and was 
first shown for the planar case for m = 1 (Jog 1999) and 
m = 2 and 3 (Jog 2000). Thus the negative disc response is 
a general feature applicable for any self-gravitating system 
subjected to an external perturbation. 

We next define a dimensionless response potential rj m 
for an azimuthal component m, at z~0 as: 



Vn 



(27) 

The negative sign is valid for odd values of m— 1,3,5 etc; 
while for even values of m = 2, 4, 6 etc the integral has a 
positive sign. The integral over /3 can be recast in terms of 
the more familiar elliptical integrals. On doing this, and sub- 
stituting back in eq.[25], we get the net response potential 
for odd values of m to be as follows: 

^°tsp(R, 4>, z , m) = —4Gcos(m<j>) 

[°° f°° p r (R',z')R'dR'dz' r _ . 

(28) 

where 



(33) 



where the response potential *P™ sp is given respectively by 
equations (28) and (32) for odd and even m values, and ^Pi is 
the perturbation potential (eq. 10), all defined at the mid- 
plane, z—0. Note that both the numerator (via the term 
H(R,z')) and denominator are proportional to the term 
GMp/D 2 in the perturbing potential. Hence the ratio r\ m 
is independent of the strength of the linear perturbation po- 
tential. 



2.3 Self-consistent disc response 

So far we have treated the disc response to the imposed 
perturbation potential alone. However the disc gravity will 
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also play a role in determining the net disc response. We next 
obtain the net, self-consistent disc response, following the 
approach of Jog (1999, 2000). A particle in the disc will be 
affected by both the imposed potential, (ip)m, and also the 
potential corresponding to the disc response to it. For a self- 
consistent case, the net potential, (ip n et)m, for the azimuthal 
number m can be written as: 



(*„et)m = (*l)m + (* 



(34) 



where ($ r esp)' m is the disc self-gravitational potential which 
corresponds to the disc response to the net potential, for 
the Fourier component m. In analogy with the direct disc 
response to the imposed potential alone (eq.[33]), the above 
can be written as: 



(*re 



= —T)„ 



(35) 



Substituting this in the previous equation, we get: 



(*„e«)m = 7-7-^— = (*lM„ 
1 + Vm 



(36) 



where <5 m (< 1) is defined to be the reduction factor for the 
m th Fourier component, and denotes the fraction by which 
the magnitude of the imposed external potential is reduced 
due to the self-consistent negative disc response. Thus, 



1 



1 + J?r7 



(37) 



Note that 5 m is always less than or equal to one by 
virtue of its definition. It tells us how strongly the disc reacts 
to counteract the effect of an external perturbation. Lower 
values of 5 denote a higher resistance due to self-gravity to 
the external perturbation. In turn, the regions with high 8 
indicate that the self gravity is weaker there and the disc is 
more susceptible to influences from the external field. From 
eq. (37), we can see that the limiting value of <5 = 1 corre- 
sponding to 77 = denotes that the disc-gravity plays no role 
in this case and it can be taken to be exposed directly to the 
external potential. On the other hand, 8 < 1, indicates a re- 
duction in the net potential due to disc self-gravity. The net 
potential determines the perturbed motion that is actually 
seen, which results in a warp for the case m = 1. 

Although the negative disc response may seem some- 
what surprising at a first glance, in reality it is a general 
feature, and we expect this to be seen in any gravitating 
system that is perturbed by an external field. It shows that 
the core of the self-gravitating system is left undisturbed 
while the outer regions suffer the consequences of the exter- 
nal perturbation potential. 



3 RESULTS 

We study the radial dependence of the resulting response 
potential for the different Fourier components m, for various 
vertical density distribution profiles in the disc. We obtain 
the plots of the reduction factor vs. radius and see where the 
minimum lies. This is next argued to denote the location of 
the onset of that particular vertical distortion in the disc. 

The perturbed motion is still described by the equation 
of motion (see eqs. [12], [13]), except that now the pertur- 
bation potential is replaced by the net potential (eq.[36]) 
that takes account of the self-consistent disc response. Thus 



the reduction factor has a minimum at say, K m in, while %f>\ 
the imposed perturbation (eq.[10]) increases linearly with 
radius. Hence the net potential has a minimum around this 
radius, and the net vertical disturbance will be seen mainly 
in the outer parts beyond this radius. Hence, we define R m i n 
to be the radius of onset of vertical distortion for the case 
m, with m = 1 for warps, as done in Saha & Jog (2006). 
The warp amplitude thus increases with radius beyond the 
onset radius. This notion of the onset radius is particularly 
applicable if the reduction factor has a sharp minimum. 

The magnitude of 8 at R m i„ and beyond is also an im- 
portant result since it tells us how strongly the disc responds 
to the external perturbation at a given radius, in view of the 
reduction due to the disc self gravity. 



3.1 Input parameters, and numerical solution 

The input parameters for the disc are: pdo, the disc central 
density; Rn, the disc scale-length; and zo, the vertical scale- 
length for the sech 2 ' n function used. The vertical density 
profile at any radius is denoted by the index n (see eq. [2]), 
where n = 1, 2 and 00 correspond to a sech 2 , sech, and an 
exponential vertical profile respectively. The halo parame- 
ters are: pno, Rc and q, or the central density, the core radius 
and the flattening parameter respectively. 

We consider a typical spiral galaxy with disc and halo 
parameters as for our Galaxy, from the model by Mera et 
al. (1998). Thus, the central surface density is 640 Mq pc~ 2 
(see Narayan & Jog 2002) , and the vertical scaleheight for an 
exponential distribution is 300 pc, hence we take the central 
mid-plane density p d0 to be = 1 M pc -3 . The disc scale- 
length R_d is taken to be 3 kpc, since this also agrees with 
that for a typical giant spiral (Binney & Tremaine 1987). 
Thus, the ratio of vertical scale-height and exponential disc 
scale-length is taken to be a constant at Zo/Rd = 0.10. Later 
these values are varied to consider the dependence on these. 

From the Mera et al. model (1998), the halo parameters 
are taken to be: puo = 0.035 M Q pc , and R c = 5 kpc. The 
halo is taken to be nearly spherical with a small oblateness, 
with the flattening parameter q = 0.95 (see Section 2.1). 
Although the dark matter halo is included for the sake of 
completeness, it has a negligible effect on the disc response, 
since the halo is less concentrated near the mid-plane than 
the disc is (eq.[18]). This is true even for a flattened halo 
as seen in cosmological simulations, see the Appendix A for 
details. 

The perturber is taken to be a satellite like the Large 
Magellanic Cloud(LMC) situated at a distance D — 48.5fepc, 
and at an angle i = 90° for simplicity. The parameter 
M p , the mass of the perturber, is taken to be 2A1O 1O M0 
(eq.[10]). In the current work, we also systematically calcu- 
late the values of the reduction factor, 8 m , for m = 1 to 10 
for a sech vertical distribution. 

Equations (28) and (32) are solved numerically using 
adaptive quadrature method for radial points spaced at 0.1 
kpc. We go upto a radial distance of 60 kpc in our calcula- 
tions for the sake of completeness. This value is beyond the 
observed extent of usual galactic discs is known to be. 
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Table 1. Warp onset radius for different vertical profiles 

m ill 



Figure X. Plot of Si, the reduction factor for m = 1 vs. radius 
R/Rd, for different vertical density profiles (sech 2 , sech and an 
exponential) in the galactic disc. As the disc profile becomes more 
flat towards the mid-plane as in the case of sech and sech 2 , the 
disc self- gravity is more important. Thus the disc resists distor- 
tion resulting in a lower reduction factor. In contrast, for the 
exponential case, the reduction factor is higher, hence the disc is 
more readily responsive to the external potential. 



3.2 Dependence on vertical density distribution 

Figure 1 contains a plot of the resulting reduction factor, 
Si for m = 1 for the various typical disc vertical density 
distributions considered (from eq.[2]). The plot shows that 
the minimum in the reduction factor and its location de- 
pend strongly on the choice of the vertical density distribu- 
tion function. The explanation for this lies in the difference 
that the exponential and sech 2 /" distributions produce. The 
sech 2 profile shows the lowest value of the minimum of the 
reduction factor, indicating the strongest negative disc re- 
sponse due to self-gravity. This is due to a higher mass con- 
centration closer to the mid-plane. In this case the distribu- 
tion is much flatter and smoothly continuous near the mid- 
plane, thereby making it more massive close to the mid-plane 
and hence more resistant to being influenced by an external 
disturbance. On the other hand, an exponential distribution 
is characterized by a sharp cusp near the mid-plane denot- 
ing less amount of matter and hence lower self-gravity. Hence 
the disc is more susceptible to the external disturbance. 

Another noteworthy result from this plot is that the 
location of the minimum depends on the vertical profile, see 
Table 1. The minimum of the reduction factor occurs at 1.9 
Rd and 2.8 Rd for a sech 2 and a sech vertical distribution 
respectively. Thus, the disc starts to respond to the external 
potential from a smaller radius which is from within the 
optical radius, hence the stellar warps will occur. Recall that 
the optical or stellar radius of a typical galaxy is observed 
to be ~ 4-5 exponential disc scalelengths, beyond which the 
intensity decreases sharply (van der Kruit & Searle 1982). In 
contrast, note the much larger radius for the onset of warps 
of 5.1 Rd for an exponential disc, close to the outer edge of 
the optical disc, as also seen in Saha & Jog (2006). Thus an 
exponential disc will mainly allow the warps to be seen only 
in HI which typically extends beyond the optical disc. 

This difference follows from the mathematical form of 
the response function and the fact that the disc gravity is 



Density profile, p(z) 
p( Z )(from eq.[2]) 

sech 2 
sech 

exponential 



1.9 
2.8 
5.1 



0.56 
0.70 
0.85 



mostly from matter close to the mid-plane for the sech 2 
case. Hence the net contribution to the response potential 
(eq. [25])is a maximum at lower radius. The minimum value 
of the reduction factor is higher for sech and the minimum 
is broader, this trend is even stronger for the exponential 
case. For the exponential case, the negative disc response 
due to the self-gravity of the disc is due to matter at larger 
z values and hence is less strongly dependent on the radial 
mass distribution. Hence the minimum in the plot of the 
reduction factor versus the radius is not pronounced, and 
the definition of Rmin is thus not so sharply defined. 

Fig. 1 thus shows that for a given perturbation poten- 
tial, a stellar disc with a steeper vertical disc density distri- 
bution has a larger onset radius for warps, and a stronger 
warp amplitude (due to a lower disc resistance as denoted 
by a higher Si value). Thus the stronger warps will have a 
steeper rise in the amplitude. This trend naturally explains 
the puzzling observation, namely "the farther away a warp 
starts, the steeper it rises", which was noted by Sanchez- 
Saavedra et al. (2003). 



3.3 Dependence on Rd and zq: 

For the sech vertical distribution, we next plot the reduction 
factor Si vs. radius for different values of the disc scalelength 
Rd for a constant vale of the disc scaleheight zo = 0.3 kpc 
(Fig. 2). Similarly, in Fig. 3 we plot the reduction factor Si 
vs. radius for different choice of zo for a constant vale of 
the disc scalelength Rd = 3 kpc. For low Rd (Fig. 2) and 
similarly for low zo (Fig. 3) the self-gravity is lower hence 
the disc follows the external perturbation well (higher Si) 
and it starts at a small radius way inside the disc. 

This dependence is further well-brought out in a sys- 
tematic study where the values for the input parameters are 
scanned, and the results obtained by solving eq.(28) in each 
case are given as contour diagrams, in the next two figures. 
We vary the values for the disc scalelength, Rd from 1 to 5 
kpc and vary the vertical scaleheight zo from 0.1 to 1 kpc, 
while keeping the central density constant at 1 Mqpc~ 3 as 
before. This then scans the behaviour of galaxies of the same 
central density but with varying mass - note that for a given 
central density, increasing Rd and zo yields a more massive 
galaxy. The resulting values of R m i„, the radius where the 
reduction factor is the minimum for m=l, vs. radius are 
plotted in Figs. 4 a and 4 b for a vertical density distribu- 
tion of type sech (top panel) and exponential (lower panel) 
respectively. In each case, the higher disc mass as given by 
higher Rd and zo leads to a higher value of R m i n in the top 
r.h.s. corner of the plot, hence the warps set in at a higher 
radius. In contrast, for a lower mass galaxy (for the lower 
l.h.s. corner in the plot), the warps set in earlier. The re- 
sults from this plot are general and can be applied to any 
galaxy. Thus, for a typical giant spiral of a given central 
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Figure 2. Plot of <5i, the reduction factor, for m = 1 vs. radius, 
R/ Rjj, shown as a function of the scalelcngth Ro for a given 
disc scaleheight zq = 0.3 kpc, for a vertical sech density profile. 
The reduction factor is lower for a higher R.£> value since the disc 
is then more massive and hence shows a higher resistance to the 
external perturbation potential. 




2 4 6 8 10 12 
R/R 

d 

Figure 3. Plot of <5i, the reduction factor for m=l vs. the radius, 
R/Rrj, shown as a function of the scaleheight zo for a given disc 
scalelength Rjy = 3 kpc, for a vertical sech density profile. The 
reduction factor is lower for a higher zo value since the disc is 
then more massive and hence shows a higher resistance to the 
external perturbation potential. 

density, if its disc scalelength and the vertical scaleheight 
are known, then this figure gives the value of the minimum 
radius beyond which warps can be seen. 

The two cases (Figs 4a and 4b) show a striking differ- 
ence, namely that for the sech profile for the vertical density 
distribution, R m i„ varies between 1.9 Rjj to 3.8 Rr>. Thus 
for all reasonable galactic parameters, the warps start to de- 
velop from within the optical disc. Hence the galaxies with 
a sech profile or with sech ' 5 as observed, will all show stel- 
lar warps. In contrast, for the steeper, exponential vertical 
density distribution (lower panel), the R m i„ covers a range 



R . /R J for expl case 

mm d 1 




Figure 4. Contour plot of R m i n the radius at which the reduction 
factor is the minimum, shown as a function of the disc scalelength 
Rb, and the scaleheight zo; for the sech vertical density profile 
(top panel) and the exponential vertical density profile (lower 
panel). The value is < 4 for most parameters for the sech case, 
implying that stellar warps would always be seen in this case. 

of larger values, from 3.5 to 7.0 Rr>. Thus over most of the 
parameter range, a disc with an exponential vertical pro- 
file cannot support stellar warps. This confirms the result in 
Saha & Jog (2006). Thus a seemingly small change in the 
vertical density profile (a sech vs. an exponential profile) 
has a drastic effect on the resulting radial onset of warps, as 
already seen in Fig. 1. 

Similarly, the results for the minimum of the reduction 
factor for m = 1 are plotted in Figs. 5a and 5b for a sech 
vertical density profile (top panel) and for an exponential 
vertical density profile (lower panel). In each of Figs. 5a 
and 5b, we see that for high zo and Kd values (top r.h.s 
corner) where the disc is more massive, the reduction factor 
is the smallest meaning the negative disc response is the 
highest. Here the disc has the maximum resistance to being 
distorted. Whereas at the lower l.h.s. corner, the disc is less 
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0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

z (kpc) 



Reduction Factor for expl case 
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Figure 5. Contour plot of the minimum of the reduction factor, 
5i , for m = 1 , shown as a function of the disc scalelength i?£> , and 
the scaleheight zq, for a sech vertical density profile (top panel) 
and for an exponential vertical density profile (lower panel). The 
reduction due to self-gravity is more significant for the sech case. 



massive, hence has less resistance due to self-gravity and 
hence responds more readily to the external perturbation. 
On comparing these two cases, one can see that the typical 
reduction factor has a smaller value for the sech vertical 
profile. This is due to a greater overall resistance by the 
self-gravity of the disc to being distorted. Thus the warps 
for the sech case would be weaker and start at lower radii. 
Thus we predict that the resulting stellar warps are weaker 
compared to the HI warps, this agrees with observations 
(Reshetnikov & Combes 1998). 



3.4 Dependence on the Fourier component, m: 

The basis of the calculation in Section 2 is the fact that 
any function can be represented as a sum of discrete Fourier 
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Figure 6. Plot of the reduction factor <5 m for different Fourier 
component m values vs. radius, R/ Rjj, for a sech vertical density 
distribution. Note that as m increases, the values of the reduction 
factor are higher and the K m i n for the minimum of the reduction 
factor shifts to lower radii. Thus there is less resistance due to 
disc self-gravity to vertical distortion. The disc is susceptible to 
higher m modes from deep within the optical disc. 



Table 2. Onset radius and reduction factor for vertical distortion, m 
m component R m „/Ri) (<5 m ) m i„ 



1 

2 
3 
5 

10 



2.S 
2.4 
2.1 
2.0 
1.9 



0.70 
0.77 
0.82 
0.84 
0.85 



components, where m — 1,2... correspond to the m th Fourier 
decomposition. It is but natural that when the density re- 
sponse is decomposed into the Fourier components, the 
higher order terms will also have an overall contribution to 
the final response potential. It is well-known that the higher 
order perturbation terms show signatures in the disc too. To 
see the behaviour of these, we calculate the reduction factor 
for a few select m upto 10 for the sech profile (Fig. 6). 

A prominent saddle-shaped m = 2 feature has also been 
seen for our Galaxy, in addition to the m=l feature namely, 
the usual integral-shaped warp (Levine, Blitz & Heiles 2006). 
Figure 6 gives the reduction factor, 82 versus R/Rd for this 
case. Similarly the higher order m results are also shown for 
m = 3,5, and 10, with an application in mind to the cor- 
rugations or scalloping observed in galaxies. We note that 
the graph for m = 2 lies higher compared to the graph of 
m = 1, and indeed all the higher m cases show a progres- 
sive trend (see Table 2). After the onset, all the Fourier m 
modes grow stronger with radius (Fig. 6), as is observed 
(Levine et al. 2006). The reduction factor values including 
the minimum for m = 2 are higher than in the m — 1 case. 
This could have significant implications as it serves to tell 
us that m = 2 signatures are seen relatively universally in a 
galaxy. It is interesting that observationally for our Galaxy 
the various modes are indeed seen over a large radial range, 
starting from a few disc scalelengths (Levine et al. 2006, Fig. 
13). This is explained naturally by our result in Fig. 6. 

Matthews & Uson (2008) have found a similar result 
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for IC 2233, where they show that the "corrugations" in the 
disc, which could correspond to high m > 1 signatures, are 
noticeable throughout the radial extent of the galactic disc, 
which are stronger in the outer regions (see their Fig. 1). 

The planar cases showed the same decreasing radial de- 
pendence of Rmm for higher m components (Jog 2000). Saha 
& Jog (2006) on the other hand got a higher R m i„ for a 
higher m =10, they probably got this wrong since they had 
not taken the correct account of the radial dependence while 
solving for the self-consistent disc response numerically. 



4 DISCUSSION 

1. Stellar warps and their detection: 

Observations show that stellar warps are common and 
occur in more than 50% of spiral galaxies(Sanchez-Saavedra, 
Battaner & Florido 1990, Reshetnikov & Combes 1998). 
These therefore must start within the Holmberg radius. In- 
deed this distinction though obvious is not often made in the 
literature- namely, a stellar warp by necessity must start 
within the optical radius. Recent systematic study of 325 
edge-on galaxies Ann & Park (2006) confirms this point, 
with a typical warp radius = 0.7 times the optical radius 
or ~ 3 disc scalelengths. Interestingly, this observed value 
agrees well with our typical onset radius of ~ 3 disc scale- 
lengths (Table 1). Thus, we have shown that a realistic, less- 
steep vertical density distribution of type sech 2 ^ n results in 
the onset of stellar warps within the optical radius as is ob- 
served in most galaxies. 

A similar value is seen for our Galaxy, where the stel- 
lar warp is shown to start from 3.1 disk scalelengths based 
on the COBE/DIRBE data (Drimmel & Spergel 2001), and 
2.4 disk scalelengths based on the 2MASS data (Lopez- 
Corredoira et al. 2002). 

The Spitzer observations of ten galaxies show the onset 
of warps lie within 3-6 disc scalelengths, thus many start 
from within the optical radius (Saha, de Jong & Holwerda 
2009). These authors treat an exponential vertical density 
profile for simplicity, and try to explain the small observed 
warp onset radii by assuming that the scaleheight increases 
with radius - while keeping the disc mass constant. Such 
flaring with radius is observed in some galaxies (de Grijs 
& Peletier 1997) and is expected for a multi-component, 
coupled, star-gas disc (Narayan & Jog 2002). However, the 
values of flaring they use are ad-hoc, and even this cannot 
explain the entire range of smaller values of R m i„ that are 
observed for their sample. Further, they use zo/Rd as a sin- 
gle thickness parameter but as we have shown (Fig. 4), the 
value of Rmin is not a simple function of this parameter. In- 
stead the small onset radii can be explained naturally as we 
have done, by using the observed less-steep vertical density 
distribution. 

2. Warp onset radius: Dependence on disc mass 

At high redshifts, the galaxy size is smaller as seen for 
the Hubble deep field sample (Elmegreen et al. 2005) with 
a typical disc scalelength of 1.5 kpc. Such a size variation 
with redshift is expected in the hierarchical evolution sce- 
nario (e.g., Steinmetz & Navarro 2002). The warp onset ra- 
dius for this sample is observed to be smaller ~ 1.4 disc 
scalelengths (Reshetnikov et al. 2002). This observed trend 
agrees exactly with our result (see Fig. 4)- namely that the 



smaller mass galaxies allow warps to develop from a smaller 
starting radius. 

At the opposite end, we predict that massive, nearby 
disc galaxies would have warp onset at a larger radius and 
thus are less likely to show a stellar warp. This would be 
tricky to confirm because the radial range over which opti- 
cal warps are seen is small (starting at 3 Rd and going up 
to 4-5 Rd), and observational data need to be analyzed to 
study this point. There is some evidence for this correlation: 
M31 has a massive disc with a scalelength of 5.4 kpc (Gee- 
han et al. 2006), and it has a stellar warp starting at radii 
larger than the isophote at 26.8 /i_g (Innanen et al. 1982). 
This is beyond the Holmberg radius, whereas the average 
onset radius of stellar warps is within this radius, see point 
1 above. 



3. Effect of nearby perturbers, and live halo : 

A real galaxy is likely to undergo close encounters with 
satellites less massive than the LMC, while our calculation is 
meant for a distant encounter with distance large compared 
to the galaxy size. Further, for simplicity, we have taken the 
halo to be rigid and the forcing frequency to be zero. The 
long-standing problem about the tidal origin of warps has 
been the resulting small amplitude (Hunter & Toomre 1968). 
This is overcome if the halo is live and the wake generated at 
the resonance points of the frequency of the perturber's mo- 
tion is included (Weinberg 1998, Tsuchiya 2002, Weinberg 
& Blitz 2006). 

The Sagittarius dwarf is about 10 times less massive 
than the LMC but is three times closer, so their direct tidal 
torques on the Galaxy (being proportional to the mass of 
the perturber/distance 3 to the perturber) are comparable. 
Hence the Sagittarius dwarf also cannot directly produce the 
Galactic warp that is observed. However, as Bailin (2003) 
has argued, the magnitude of the tidal field due to the wake 
generated in the halo by the Sagittarius dwarf could also ex- 
plain magnitude of the warp seen in the Galaxy. This needs 
to be checked by simulations. 

Yet another pathway to create vertical perturbations 
would be to have an even smaller mass perturber, like the 
subhalo, come even closer to ~ 5 — 10 kpc. This possibility 
has been studied by Chakrabarti & Blitz (2009) via simula- 
tions, who show that this can explain the HI amplitudes of 
warps and the higher order vertical modes as observed by 
Levine et al. (2006). 

In the present paper we have not attempted to obtain 
the actual warp amplitude. However, the concept of negative 
disc response studied here would still apply in these general 
cases. That is, due to its strong self-gravity, the inner disc 
region would resist being distorted by vertical perturbations. 

4- Dynamical implications: The present paper shows the 
surprisingly strong dependence of the resulting warp radius 
on the vertical density distribution in the disc. This is be- 
cause the disc self-gravity which decides the warp radius de- 
pends crucially on the vertical disc distribution close to the 
mid-plane. A similar strong dependence on the disc vertical 
distribution may also affect other dynamical studies such as 
the vertical heating due to tidal encounters, or the bending 
instabilities. We plan to look at these in future papers. 
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5 CONCLUSIONS 

We study the self-consistent, linear response of a galac- 
tic disc to vertical perturbations arising due to a tidal en- 
counter. We show that the vertical disc density distribution 
has a surprisingly strong effect on the radius at which disc 
starts to show the vertical distortions. In retrospect this is 
physically understandable since the disc resists distortion in 
the inner parts due to its self-gravity, which in turn depends 
crucially on the vertical mass distribution close to the mid- 
plane. A fiat-core profile like the sech as is observed results 
in a small radius for onset of warps of ~ 3 disc scalelengths. 
We can thus naturally explain why the radius of onset of 
stellar warps is within the optical disc as observed in most 
galaxies. The results for the radius for warp onset are given 
as a function of the disc scalelength and disc scaleheight, 
and are presented as contour plots, which can therefore be 
applied to any galaxy. 

We show that the higher m distortions develop from 
even smaller radii ~ 2 disc scalelengths. Hence the higher 
mode corrugations can be seen over a large radial range, in 
agreement with what is observed for the Galaxy (Levine et 
al. 2006), and for IC2233 (Matthews & Uson 2008). 

Acknowledgments: We would like to thank the anony- 
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us to include the effect of the shape of the dark matter halo 
on the warp. We also thank Lynn Matthews for useful com- 
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Appendix A: Weak dependence of warp on the shape 
of the dark matter halo 

We have included a spherical dark matter halo in the 
calculations (Sections 2 and 3) for the sake of completeness. 
We show next that it has a negligible effect on the vertical 
response of the disk including the warp onset radius. In this 
Appendix, we also investigate the possible effect of the shape 
of the dark matter halo on the warp. Cosmological N-body 
simulations show that the dark matter halos are not spher- 
ical but instead tend to be flattened, with the typical value 
of q = c/a, the vertical to planar axes ratio of 0.6 (Bailin & 
Steinmetz 2005) or 0.7 (Bett et al. 2007; Kuhlen, Diemand 
& Madau 2007). 

In this paper, we have studied the self-consistent re- 
sponse of a galactic disc to an external perturbation, and 
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Figure 7. Plot of the reduction factor 81 for disc-alone, and for 
disc-plus-halo with different oblateness of the halo , with q = 0.95 
(spherical) to q=0.4 (highly flattened halo) vs. radius, R/ Rd, 
for a sech vertical density distribution. The spherical halo has a 
negligible effect on the net disc response. Even a flattened halo 
with the axis ratio, q =0.7 as seen in cosmological simulations, 
plays a minor role in determining the disk response, changing the 
radius of onset of warps by just 10%. 



Table 3. Warp onset radius for different halo flattening 



Component 




Rmzn/R.D 




disc- alone 




2.6 


0.70 


disc-plus-halo (q= 


=0.95, spherical) 


2.7 


0.71 


disc-plus-halo (q= 


=0.7, flattened) 


2.9 


0.72 


disc- plus- halo (q= 


=0.4, very flattened) 


3.2 


0.74 



on the radius of onset of warps given by the location of the 
minimum in the reduction factor curve (see Table 3). 

The shape of the halo is also shown to play a minor 
role in determining the disc response. A typical oblate halo, 
with a flattening of q — 0.7 as seen in the cosmological 
simulations, changes the radius of onset of warps by just 
10% (Table 3). Interestingly, even a highly flattened halo 
with q = 0.4, which occurs at the most oblate end of the 
halo distribution found in simulations (e.g., Bett et al. 2007, 
Bailin & Steinmetz 2005), has a small effect on the radius of 
warps. Thus over the entire range of flattening of the halos 
seen in cosmological simulations, the halo has a minor effect 
on the determination of the warp. 



the halo affects this calculation in an indirect way. The ex- 
ternal, tidal perturbation affects both the disc and the halo. 
The self-consistent disc response is obtained by treating its 
response both to the external potential, plus the response 
potential corresponding to the disc and halo response to the 
net perturbation (eq. [34]). Thus the inclusion of halo enters 
in the calculation of the frequencies (eq. [14]), and the net 
response density (eq. [18]). Of these, the frequencies v 2 and 
Q 2 are affected by the halo, especially in the outer parts and 
for an oblate halo. On the other hand, the response density 
of the halo is small ~ 10 -3 x the disc response density, hence 
it does not affect the self-consistent disc response. This is be- 
cause the halo, even an oblate one, is much less concentrated 
towards the disc mid-plane. 

We solve eq.(28) (see Section 3.1) and obtain the reduc- 
tion factor Si vs. radius for a sech vertical profile of the disc, 
for a disc-alone case, and the disc plus halo cases for q = 0.95 
(spherical case), q=0.7 (the typical value seen in cosmolog- 
ical simulations), and q=0.4 (an extremely flattened oblate 
halo). The disc and the halo parameters for the spherical 
case are as given in Section 3.1. For the flattened halo cases 
the relation between the central density p c and the core ra- 
dius R c are obtained in terms of those for a spherical halo, 
assuming a constant halo mass, and using the relations de- 
rived in Narayan et al. (2005, see their eq. [7]). The results 
are plotted in Fig. 7. 

The most striking result is that the inclusion of a spheri- 
cal halo has a negligible effect on the disc response. The low- 
est two curves are nearly identical. The weak dependence on 
the halo may seem surprising since the halo is known to be 
dominant in the outer parts. However, there are two reasons 
for this: first, the disc mass and the halo mass are compa- 
rable upto the region of interest, namely R m i„ ; second, the 
halo is much less concentrated towards the disc mid-plane. 
Thus, the inclusion of a spherical halo has a negligible effect 



